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Abstract. The effects of quadratic order terms in the dispersion matrix near a mode 
conversion are considered. It is shown that including the corrections due to these quadratic 
' terms gives a better matching between the local solution in the mode conversion region. 
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and the far-field WKB solutions for the incoming and outgoing waves. This matching is 
demonstrated by comparison of the asymptotic solution with a numerical solution for a simple 
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Q ' one-dimensional conversion. This procedure for obtaining the corrections due to quadratic 

' , order terms can be extended to arbitrary order and, in principle, an outline for performing such 

an extension is given. 
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^ ! 1. Introduction and Motivation 

oo 

O ■ The resonant interaction of linear waves in a nonuniform, or time dependent, medium is a 
phenomenon of great interest in a wide variety of fields. There is a large literature going 
^ ■ back many years which we will not attempt to survey here (see, for example, the literature 
o3 ■ cited in uj). In this paper, because the resonance is local, either in time (for time-dependent 
backgrounds) or in space (for nonuniform backgrounds), we assume that away from the 
resonance the evolution is well described using adiabatic methods (e.g. WKB methods). 
Hence, the problem reduces to ensuring a good matching between incoming and outgoing 
WKB solutions. The resonance is dealt with using an appropriate local evolution equation 
that is simpler to solve than the full equation. Often this is done by simply linearizing the 
background dependence of the coupled evolution equations in the immediate vicinity of the 
resonance, leading to a Landau-Zener-type model. The local solutions can be expressed in 
terms of parabolic cylinder functions which are then matched to the incoming and outgoing 
WKB solutions. This level of approximation to the local dynamics has the defect that the 
parabolic cylinder functions, while capturing the local jump in wave amplitude due to the 
resonant interaction, fail to capture the slower amplitude variation away from the resonance. 
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This amplitude variation, which is due to action conservation in a nonuniform or nonstationary 
background medium, is well modeled by the WKB solutions. Hence, the asymptotic matching 
region, where the true solution is well modeled by both WKB and parabolic cylinder solutions, 
can be quite small, leading to problems for automated matching procedures. Here we consider 
how to systematically account for higher order effects in order to improve the robustness of 
the matching procedure. Including the quadratic order corrections, for example, dramatically 
enlarges the matching region. Quadratic order effects have already been considered by many 
authors. We mention, for example, Delos and Thorston |l2l, Swanson and Friedland, et 

al. The new contribution of the present paper is that we consider higher order corrections 
from a phase space perspective. This has the advantage of casting the problem into the 
simplest possible form, which is universal in character, while preserving all of the essential 
features of the particular problem. By analogy with problems in particle mechanics, the 
phase space approach allows a wide variety of changes of representation for simplifying the 
problem as much as possible and cleanly separates the inherent complexity of the problem 
from complexity which is simply due to a choice of representation. While this is worked out 
in detail only for the case of quadratic order corrections, in the Appendices we also sketch an 
argument describing how the results might be formally extended to arbitrary order. Another 
great advantage of phase space methods is that they are applicable to multidimensional 
systems, which is an area of ongoing research. 

Of central importance in the phase space theory of mode conversion is the concept of 
a normal form for the local equation. This is, in a very well-defined sense, the simplest 
representation of the local dynamics. The natural physical interpretation is that it is the 
representation that puts the 'uncoupled' wave operators on the diagonal, with the off-diagonals 
giving the coupling. This change of representation to normal form should be carried out by 
an adiabatic transformation, meaning that the polarization basis used to construct the local 
representations should be well behaved everywhere and not just away from the resonance. 
Normal forms for one-dimensional conversion, like we consider here, were studied in great 
detail by Littlejohn and Flynn In dimensions higher than one we mention the work of 
Littlejohn and Flynn fT\, Braam and Duistermaat (HIIH, Colin de Verdiere lfTOl[TT| . Kammerer 
ifTIl and Tracy and Kaufman [fT3l. 

Because the most natural arena to view WKB theory is phase space, we wish to view 
conversion as a phase space phenomenon and develop methods based upon the geometric 
ideas of Maslov theory. In previous work [[H [6l [TH [151 it was shown that phase space 
techniques can be used to solve wave problems exhibiting mode conversion between modes of 
two different polarizations. Such multicomponent problems can be written initially in x 
matrix operator form. Using the congruent reduction procedure developed by Friedland and 
Kaufman in lfT6ll . the system can be reduced to a scalar wave equation away from resonances, 
leading to traditional WKB methods. In the vicinity of mode conversion, the congruent 
reduction procedure can be used to reduce the problem to a 2 x 2 matrix form governing 
the two interacting modes. The advantage of the congruent reduction procedure, which uses 
general congruence transformations rather than unitary transformations, is that the change 
of polarization basis used to carry it out can be smooth everywhere. Diagonalization using 
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unitary operators does not have this property due to the near-denegeracy of eigenvalues, which 
is the essential characteristic of mode conversion. 

As an example of a physical problem which exhibits mode conversion, consider the 
equations for the propagation of an electromagnetic wave in a cold plasma. From Maxwell's 
equations we have 

V X V X - ^^K ■E = Q. (1) 
The dielectric tensor K for a cold plasma can be shown to be ifTTl 
/ S -iD \ 

¥.= iD S Q , (2) 
^0 P y 

where the functions S, D, and P depend on the density of the plasma, and the background 
magnetic field. Spatial variations in these plasma parameters lead to spatial variations in the 
nature of propagating solutions, and in certain conditions |[T5l , this problem will exhibit mode 
conversion. 

Returning to a generic mode conversion problem, we consider the problem at a given 
frequency, and write the two equations for the coupled wave modes together, 

D(x,-i9^;cu) ■ = 0. (3) 

While this equation works for waves in multiple spatial dimensions, in this paper we limit 
our analysis to the case where x is one-dimensional. Additionally, we will suppress the 
frequency dependence for brevity of notation. Using the symbol calculus ifTSl [T9l . we can 
define the symbol of the wave operator as a matrix-valued function on wave phase space, 
D(x, k). Here, the variable k corresponds to the operator —id^, and products of conjugate 
variables correspond to symmetrized operators (e.g., the symbol xk maps to the operator 
-i{xd^ + d-,x)/2). 

In the vicinity of a mode conversion, there are two roots of the dispersion relation 
det(D(x, k)) = 0. These two curves in phase space locally have a hyperbolic structure (an 
"avoided crossing", see Figure ©). Linearizing the x and k dependence of the dispersion 
matrix about the center of the hyperbola, and then converting this linearized symbol back 
to an operator, gives a set of coupled equations which can be solved for the local wave 
fields. Matching these local solutions onto uncoupled WKB solutions (which are a good 
approximation to the solutions far from the mode conversion region) gives transmission and 
conversion coefficients for the incoming and outgoing waves. These coefficients can be used 
to treat the mode conversion as a ray-splitting process, where the incoming ray is split into 
two outgoing rays, one for each mode. 

This ray-splitting approach captures the jump in amplitude caused by the coupling 
between the two modes at linear order in phase space variables. However, higher order terms 
in the wave equation can lead to additional effects. For example, the amplitude variation 
familiar from WKB theory is not captured by the linear order solution. This could cause 
difficulties when attempting to match the local wave fields onto the incoming and outgoing 
WKB solutions. Figure ^ illustrates this effect by comparing numerical solutions to the 
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Figure 1. The phase space structure of a typical "avoided crossing" mode conversion. The 
hyperbolic dispersion curves are the solid lines, and the dashed lines are the dispersion curves 
for the uncoupled modes. 




Figure 2. Dispersion surfaces for the uncoupled WKB modes, defined by Z^n (g, p) ~ Q and 
^22 (9, p) = 0, cross at the mode conversion point. The coupled dispersion surface, defined 
by solving det(D(q,p)) = 0, has a hyperbolic structure in the mode conversion region. Its 
branches asymptote to the uncoupled dispersion surfaces. Quadratic terms in the expansion 
of Dii{q,p) and D22{q,p) about the mode conversion point cause the dispersion surfaces to 
curve. The coordinates are described in Section lTT] 
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linear order solutions. In this example, the local solution captures the jump in amplitude at 
the mode conversion, but misses the slow WKB amplitude variation. This limits the matching 
region to a small range right near the mode conversion, which could make numerical ray- 
tracing algorithms somewhat unstable. However, this example also suggests that it may be 
possible to calculate the local solution and include the effects of the higher order terms. As 
will be shown, one of these effects is an amplitude variation which matches the slow WKB 
amplitude variation. 

In this paper we will first briefly review the linearization procedure which leads to the 
parabolic cylinder functions as local solutions in the mode conversion region. We will then 
consider the effect of adding generic quadratic terms to the wave equation. These terms 
will modify the uncoupled dispersion relations (which are given by the diagonal elements of 
the dispersion matrix, after it has been transformed into a "normal form"), which changes 
the far- field WKB solutions for the incoming and outgoing waves. When the new quadratic 
terms are added to the coupled equations, they will give us a new local solution for the wave 
fields. This local solution will contain both non-propagating "near- field" contributions and 
propagating contributions to the original, linear order, local solution. The near-field terms do 
not propagate, and therefore do not affect the matching. The other modifications, however, 
are phase and amplitude corrections which make the local solution better match the far-field 
WKB solutions. We find that the quadratic order terms do not modify the S'-matrix (WKB 
connection) coefficients. Last, we give a comparison with numerical solutions for a simple 
example. 

There are three appendices to this paper. In Appendix A we show how to transform the 
dispersion matrix for the one-dimensional problem into normal form, through second order 
in phase space variables. In Appendix B, we argue that the calculations of this paper can 
be extended to higher order in the phase space variables, since it is possible to transform 
the dispersion matrix into normal form at any order. Finally, in AppendixCj we calculate 
corrections to the normal form transformation which are due to the Moyal star product. We 
show that, while these effects can be calculated, they introduce a higher order correction, and 
therefore can be neglected at the quadratic order we consider in this paper. 



2. Qualitative Discussion of Results 

There are three primary results presented in this paper. These are 

(i) New local solutions are found, which include the effects of quadratic order terms in the 
dispersion matrix 

(ii) The matching region is expanded, 

(iii) The transmission the conversion coefficients are unchanged through 0{e) where e is a 
formal small parameter associated with the corrections. 

For the first result, we Taylor expand the dispersion matrix out to second order in 
phase space variables, and truncate. We then choose a particular representation, and use 
it to convert the truncated dispersion matrix into a pair of coupled equations for the fields. 
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Figure 3. Comparison of WKB (black dashed), numerical (gray), and linear order local (red) 
solutions. 



This approximate wave equation is valid near the mode conversion, and solutions to this 
equation are new local solutions. The new local solutions have the form of the local solution 
obtained from truncation of the dispersion matrix at linear order, but with corrections due to 
the quadratic order terms. 

The expansion of the matching region can be seen most clearly in the comparison of the 
new local fields to the WKB approximations which are valid far from the mode conversion. In 
Figure ([3]) the linear order local solutions are plotted along with the WKB approximations and 
"exact" numerical solutions. The local solutions can be used for matching, but higher order 
terms in the equations result in amplitude and phase variations in the WKB solutions that 
are not in the local solutions. This means the matching must be done close the the mode 
conversion point in order to obtain reasonable results. In Figure (U), the corrected local 
solutions are plotted instead of the linear order solutions. The quadratic order corrections 
clearly make the new local solutions match the numerical and WKB solutions over a much 
larger range. This means that the region in which the matching can be accurately performed 
has now been greatly expanded. 

That the transmission and conversion coefficients are unchanged by the second order 
terms will be seen in Section 14. 3[ where the functional forms of the new local solutions and 
the WKB solutions are matched asymptotically. The implication of this result is that the 
previously obtained transmission and conversion coefficients can still be used in analysis or 
numerical routines, without modifications due to quadratic order terms. 
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Figure 4. Comparison of WKB (black dashed), numerical (gray), and quadratic order local 
(red) solutions. Since the numerical solutions are almost entirely covered by the local 
solutions, the subfigures show the difference between the numerical and the quadratic order 
local solutions. The parameters for this simulation are e = 1, ai = 3 x 10^'^, bi = 7.44 x 10"'^, 
ci 1 X 10"'', 02 = 3 X 10-3, 62 = 7.44 x lO'^ and C2 = -5 x 10"*. 



A further result which will be discussed in Appendix B is that the ideas and analysis used 
in this paper can be used to obtain further corrections to the local fields, order by order in the 
phase space variables. The Taylor's series for the dispersion matrix can be can be truncated 
at some order, and then thrown into normal form. This new approximation to the dispersion 
matrix can then be used to write equations for the local fields. Solutions to these equations 
will give new expressions for the local fields, which will contain corrections out to the order 
where the series was truncated. 



3. Linearization of the coupled system 

This review of the solution of the mode conversion problem follows the solution described 
in Reference Q. We give this review of the solution in some detail because the solution to 
the quadratic order problem proceeds along the same lines. The first step of the solution is to 
linearize the symbol of the wave operator about the mode conversion point. Then, transform 
the linearized symbol via a change of polarization basis and a linear canonical transformation 
of phase space, in order to simplify the symbol as much as possible. Convert this approximate 
symbol back into an operator, which gives a new "local" wave equation which can be solved 




Figure 5. The dispersion curves in three different coordinate systems: (a) the physical 
coordinates {x, k), (b) the symmetrized coordinates {Q, P) of Flynn and Littlejohn [61, and 
(c) the coordinates {q,p) of Tracy and Kaufman ifTSl . 



analytically. The solutions can then be analyzed in various representations, which correspond 
to different choices of coordinates in phase space. The transmission and reflection coefficients 
are then derived by matching the local solution onto the far-field WKB solutions. 

3.1. Coordinates Used in this Calculation 

By performing this calculation in phase space, the underlying symmetry between 
configuration coordinates and wave numbers becomes evident, and the techniques of linear 
canonical transformations can be used to simplify the problem. This can lead to notational 
difficulties, however, since any particular choice of linear canonical transformation of the 
phase space variables can lead to a new notation for the variable's names. In this paper, 
there will be primarily three coordinate systems used to label the phase space variables, as 
illustrated in Figure ©. The first set of coordinates, (x, k), is used for the physical coordinates 
in which the problem is naturally written. In these coordinates, the hyperbolic structure of 
the dispersion surfaces near a mode conversion will be oriented at some arbitrary angle in 
phase space [see Figure ([5^)]. The second set of coordinates, denoted {Q,P), is the pair 
of symmetrized coordinates used by Flynn and Littlejohn in BH. In these coordinates the 
asymptotes of the hyperbolic avoided crossing are oriented on the diagonals [see Figure dSh)]. 
Also, these variables have been normalized so that the diagonal elements of the linearized 
dispersion matrix form a canonical pair. The final set of coordinates is that used by Tracy 
and Kaufman in [|18|. In these coordinates, denoted {q,p), the axes are the asymptotes of the 
avoided crossing [see Figure dSt)]. The diagonal elements of the linearized dispersion matrix 
are again normalized to form a canonical pair. An advantage of this choice of coordinates 
is that it puts the linearized dispersion matrix into a kind of normal form, where the two 
diagonal elements take the simple form Dii(q,p) = —p and D22{q,p) = q- As we will see, 
the equations for the local fields are particularly simple to solve in these coordinates. 

Associated with each of these three different choices of coordinates is a different 
representation of the fields. The linear canonical transformation from one pair of coordinates 
to another induces a metaplectic transformation of the fields (which can be thought of as a 
generalization of the Fourier transform). For example, the linear canonical transformation 
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from the variables (g, p) to {Q, P) is 

The associated metaplectic transformation is 

^(g) = / e*^^(Q'^)^(g)c/g, (5) 




where Fi{Q, q) is the mixed-variable generating function for the canonical transformation: 

F^{Q,q) = ]^{Q^ - 2V2Qq + q^). (6) 

These equations illustrate the relationship of the metaplectic transformation to the Fourier 
transform. If the transformation matrix in Equation (H)) is squared, it would give a rotation by 
7r/2 in phase space, and the metaplectic transformation in Equation ([5]) (iterated twice ) would 
become a Fourier transformation. So, Equation ^ can be thought of as the "square root" of a 
Fourier transformation. 

3.2. Linearization and Solution 

The Taylor's series expansion of the symbol of D about the mode conversion point (x*, k^) is 
given by 

D(x, k) = D(x„ k,) + {x - X,)— + {k - k^)— + ^ Y^^~d]^ "' 

where all derivatives are evaluated at the conversion point (x*, k^). First shift the origin 
in phase space to (x*, k^). The shift of origin in x is performed by the change of variables 

rp V /y» rp • 

Jb ~ Jii JU ^ J 

i)\x)=i){x-x,). (8) 
The shift in k is accomplished by multiplication by a phase: 

^"(x) = e'^'''^{x). (9) 
That this is a shift in k can be seen by computing the Fourier transform of il)"{x): 

[J^^}j"]{k) = j e-'^''i)"{x)dx = j e-'^''+'^*''i){x)dx=[Ti)]{k-K). (10) 

Now truncate the series in Equation (|7]) at linear order in phase space variables. As 
shown in Reference the resulting matrix can be simplified by choosing the polarization 
basis which makes its off-diagonal elements constant: 

D(..t)=f ^"(f*> „,",,). (11) 

We now have the situation shown in Figure dS^). The dispersion surface (iet(D) (x, k) = has 
a locally hyperbolic structure, and the center of the hyperbola is at the origin in phase space. 



Quadratic corrections to the metaplectic formulation of resonant mode conversion 



10 



The asymptotes of the hyperbola are given by the zeros of the diagonal elements of D(x, k). 
These diagonal elements can be simplified by performing a linear canonical transformation 
of the phase space coordinates {x,k) to the coordinates of Figure dSj;). This puts 

the asymptotes of the hyperbola on the q and p axes. With this choice of linear canonical 
transformation, the linearized dispersion matrix becomes 

D(g,p)= ( :f ], (12) 
where the coupling f] includes the Poisson bracket of the diagonals as a normalization; 

^ = TD—rrT\ • ^^^^ 

The complex constant f] is the coupling between the two modes. In the limit f] ^ 0, the 
two modes would become uncoupled. The dispersion relations for the two coupled modes 
are given by setting the eigenvalues of the dispersion matrix equal to zero. For f/ — 0, the 
polarization basis chosen in Equation (fT2l) diagonalizes the dispersion matrix, and therefore 
the diagonal elements are the dispersion functions for the uncoupled modes. For nonzero fj 
the diagonals are no longer equal to the eigenvalues. However, far from the mode conversion 
region the diagonals are approximately equal to the eigenvalues, and so provide a good 
approximation to the dispersion functions for the two modes. This means that the WKB 
method can be used to construct approximate solutions, using either the diagonal elements of 
the dispersion matrix (in the "normal form" given by Equation (fT2l) ) or the eigenvalues of the 
dispersion matrix. We refer to these two types of approximate solutions as the uncoupled and 
coupled WKB solutions, respectively. 

The approximate dispersion matrix in Equation (fT2l) can now be converted back from 
a symbol into an operator, and we have an equation for the local wave field in the q 
representation, 

= 0. (14) 




V 

Notice that the second row of this matrix equation is now an algebraic equation instead of a 
differential equation. Solving this algebraic equation for ijj2{q) in terms of ipi{q) and inserting 
the resulting expression into the first row gives a first-order differential equation for ^i(g). 
Integrating this equation, we obtain the following solution: 

= I _A^^-i\fi\Hng 1 = 1 _^f^*q-i\v\^-l I ' (15) 

where A is a constant of integration. This expression must be evaluated with care for negative 
values of q because of the branch in the complex logarithm. If we take the branch to be just 
below the negative axis, then ln(g) = ln(|g|) + in and the solution can be written in terms of 
the absolute value of q: 

ni) = { i = i h a^) 
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where the function /(g) is defined for real q as 

{I It for g < 
1 for g > 

in order to deal with the singularity at g = 0. The variable r = e"'^!'^'^ is the transmission 
coefficient. If we absorb the factor of 1/r into the definition of A, then we can see that the 
absolute value of decreases by a factor of r as we go from negative g to positive g. For 
this solution, the energy lost from the upper channel is converted into energy outgoing in the 
lower channel. 

The solution given in Equations (fT6l) and (flTI) is only one of two possible types of 
solutions for the local fields. The other solution can be found by placing the branch cut just 
above the negative g axis, instead of below it. In this case, ln(g) = ln(|g|) —in for g < 0. This 
means that now the amplitude of ipf^ increases by a factor of 1 /r as we go from negative to 
positive g. This solution corresponds to the case where energy is coming in on both the upper 
and lower channel, but the amplitudes and phases are just right so that all of the energy leaves 
in the upper channel. 



3.3. Alternative Representations 

The form of the solution given in Equation (fT6l) is particularly convenient for solving the 
linearized system of equations, and calculating the transmission and conversion coefficients. 
However, the coordinates (g,p) of Equation (fT6l) are related to the physical coordinates 
(x, k) by some linear canonical transformation. Therefore, to find the solution in the x 
representation requires computing a metaplectic transformation of the solution \E'(g). An 
additional inconvenience of using the g representation is that the solution for the lower channel 
has a singularity at g = 0. While this is not too surprising since the dispersion manifold for 
this mode is the p axis, it makes analysis of this function tricky. 

As an example of how to transform the solution in Equation (fT6l) into a different 
representation, we will transform to the variables (Q, P) described in Section 13.11 This 
transformation induces the metaplectic transformation given in Equation ([5]). A table of 
integrals or a computer algebra system like Maple can be used to evaluate this integral, 
resulting in the Q representation of the upper channel being written in terms of the parabolic 
cylinder function U{a,Q): 

^;o)(Q) = j e^^i(Q''')^i°)(g) dq (18) 

= Ae^""'/^ U - 1/2, -(1 + i)Q) . (19) 

Here, A is a complex amplitude, whose value is set when we match this local function to 
the incoming wave. It is related to the amplitude A of Equation (fT6l) by a constant factor; 
A = -i^e'^^^^^I^A. 

Because of the singular nature of 4'2^\q), the metaplectic integrals to convert this to a 
different representation become difficult to evaluate. In [HI, Tracy et al. compute the Fourier 
transform of iP^2\q) using contour integrals, which gives the p representation of the lower 
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channel. An alternative approach is to use the Q representation of the wave equation to write 
il)f\Q) in terms of the parabolic cylinder function and its derivatives. 

^f\Q) = J e^^^^^'"^ (^^r(g)) dq (20) 

= 1 f {tdge'^'^^''^)^f{q)dq (21) 

= (^1=(Q- z9q) e^^^^Q'")^ 4'\q) dq (22) 

= ^{Q - zdQ)ijf'\Q) (23) 

= -Afi*U {i\fi\^ + l/2,-{l + i)Q) . (24) 



Recurrence relations for the parabolic cylinder function provide the simplification to the 
last line above. This representation of the solution can be easily compared to numerical 
simulations of the original system of equations, since methods for calculating the parabolic 
cylinder functions are readily available. See Figure dS]). 



3.4. Transmission and Conversion Coefficients 



The transmission and conversion coefficients for this problem are derived from asymptotically 
matching the local solution in Equation (fT6l) to the far-field WKB solutions, which are defined 
using the eigenvalues of the dispersion matrix as the ray hamiltonians. Because of the rapid 
variation of the eigenvectors in the mode conversion region, the WKB solution breaks down 
there. However, there is a region where both the local and WKB solutions are valid, and 
have the same functional form. By matching the solutions in this region, a globally valid 
approximate solution can be formed. 

The eigenvalues of the dispersion matrix in Equation (fT2l) are 



A- 



q-p 



± 



q + p 



(25) 



2 V V 2 

The choice of which eigenvalue to consider is determined by the sign of q. For g < 0, there is 
a solution to A+ = 0, while for g > 0, a zero of the other eigenvalue A„ exists [see Fig. Qc))]. 
Setting A = and solving for p{q) gives the WKB phase: 

0(9)= / Piq')dq= [ (-^^ dq' = -\fi\'^\n{q). 



The WKB amplitude is found by taking the p derivative of the eigenvalue: 



A{q) 



dp 



-1/2 



p=^\fi\2/q 

Evaluation of the p derivative, and simplification of the result, gives 



A{q) 



O 



I ~I4 



-1/2 



2^ 



O 



I ~I4 



(26) 
(27) 

(28) 
(29) 
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where we have used the fact that we are interested in the eigenvalue = \-sgn{q)- Together, 
the WKB amplitude and phase give the WKB solution for the upper channel: 



WKB/ N _ „-i|77|2ln{9) 



2g2 ^ V 9^ 



(30) 



The transmission coefficient is obtained by matching this WKB solution to the local solution. 
However, far from the mode conversion, where the WKB approximation is valid and the 
matching can be performed, we have ^ |f/p. This means that the C (^^) term can be 
dropped from the amplitude. We are then left with an expression that has the same functional 
form as the local solution given in Equation (fT6l) . This confirms our previous identification of 
r = e^^l''!' as the transmission coefficient. 

The derivation of the conversion coefficient is complicated by the need to evaluate the 
lower channel in the p representation. The integral in the Fourier transform which changes 

(q) into the p representation must be treated as a contour integral because of the singularity 
at g = 0. The details of this calculation are given in Reference |[T1, and the result is that (p) 
is zero for negative p, while for positive p it is 



4°'(P) = --^p'""'- <31) 

Comparing this to t/jf'^ (q) for negative q we find that the conversion coefficient is 



P- - (32) 



The incoming and outgoing amplitudes are therefore related by a scattering matrix: 



where the amplitude functions are defined as in Reference IfTSi 




(33) 



^f\±q) = \qi\'<\'ijf\±q) (34) 
(±p) = b|-l'^l'4°H±p). (35) 



4. Extension to Higher Order 

In general, the Taylor's series expansion of the dispersion matrix in Equation (|7]) will contain 
terms of higher order than the linear terms kept in the approximation of the previous section. 
For example, second order derivatives, such as the curl of the curl in Equation ([U), will lead to 
terms quadratic in p in the dispersion matrix. Such higher order terms will have an effect on 
both the far-field WKB solutions, and on the local solution. If we can calculate these effects, 
then we can obtain a better match of the local solution to the incoming and outgoing WKB 
solutions. This will allow us calculate any corrections that there may be to the scattering 
coefficients. 
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As shown in [Appendix A[ the quadratic approximation to the dispersion matrix can be 
put into the normal form 

The diagonal elements of this normal form matrix include arbitrary quadratic terms in addition 
to the linear terms which were described previously, while the off-diagonal terms are constant 
through quadratic order in the phase space variables. Introducing the small parameter e to 
keep track of the quadratic terms we obtain 

Dn{q,p) = -p + e{aip^ + hipq + Ciq^) (37) 

and 

D22{q,p) = q + e(c2P^ + hpq + 025^)- (38) 

This second order matrix valued function on phase space can now be converted into a 
pair of coupled second order differential equations for the two modes, in the q representation: 

= 0, (39) 




where the operators on the diagonals are 

Dn = id, + e{a^{-id,f - '^{d.q + qO,) + Cig^), (40) 

and 

D22 = q + e{c2{-id,f - '^{dgq + qd,) + a^q^). (41) 

As seen in the previous section, the lowest order solution for the second channel has a 
singularity at the origin when it is written in the g-representation, and therefore was Fourier 
transformed to obtain the p-representation given in Equation (|3T| ). Anticipating that this form 
of solution will persist even in the presence of the new (9(e) terms, we will analyze the lower 
channel in Fourier space. The Fourier transform of these equations gives the p-representation 
of the equations: 

D'li V \ ( ^i(p) 



fj* D',, J \ Mp) 
which has diagonal elements 



0, (42) 



D'n = -p + e{aip^ + ^ W + ^pP) + (43) 



and 



= tdp + eia^iidpf + '-^{pdp + dpp) + C2/). (44) 

We will proceed by first finding the new far-field WKB solutions, and then finding 
the coupled local solution. Then we will examine the local solution and match it onto the 
WKB solutions. Doing so will show that no new corrections to the scattering coefficients 
are necessary. Finally, we compare the new local solutions to numerical solutions in order to 
show the improved matching. 
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d) 

















Figure 6. Dispersion curves for the uncoupled modes, showing the effect of the higher order 
terms (e = 1/5). The red curves are solutions to Dn = 0, and the blue curves are solutions 
to £'22 = 0. The original crossing is at q = p = 0. For all plots 02 = ^2 = C2 = 0. The 
remaining parameter values used for these plots are: (a) ai = 1, bi = c\ =0: This term 
introduces a new branch to the dispersion surface D\\ = at p = l/oie. As e 0, this 
branch moves off to infinity, (b) ci ~ 1, ai = 61 = 0: This term introduces a curvature into 
the D\\ dispersion surface. In the p representation, an uncoupled solution in the upper channel 
would now look like an Airy function, (c) 61 = 1, ai = ci = 0: Although not obvious from 
the dispersion curves, this term introduces an amplitude variation which has a square-root 
form, and which matches the WKB amplitude variation due to action conservation. This term 
also introduces a resonance in the upper channel at g = l/6ie. (d) 61 = 1, ci = 0.1, ai = 0: 
In general, there will be several higher order terms, all of which effect the dispersion surfaces, 
leading to phase shifts, amplitude variations, and resonances. 



4.1. Effects of Quadratic terms on the WKB solutions 

The calculation of the WKB solution proceeds as in Section I3.4[ where the WKB solution 
was calculated to linear order in the phase space variables. Introduction of the higher order 
terms will modify the eigenvalues of the dispersion matrix, and introduce corrections of (9(e) 
into the WKB solutions. 

First, we can use the equation detD(g,p) = to find the WKB phase, since the 
determinant is only zero when one of the eigenvalues is zero. When we include the second 
order terms in the dispersion matrix, this equation becomes 

detD(g,p) = {-p + eVi{q,p)){q + eV2iq,p)) - = 0, (45) 
where the second order terms are 

Vi{q,p) = aip^ + biqp + ciq'^, V2{q,p) = a2q^ + b2qp + C2p'^ . (46) 
Expand the solution p(g) in powers of e 

p{q) = Po{q) + m{Q) + ■ ■ ■ (47) 
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and insert this into Equation (|45] ). The terms of order e° give 

0, 



pq — |?7|^ 



16 



(48) 



and so the first term in p{q) is 
Po(g) = 



?7p 



The order e terms in Equation (l45l) can be solved for pi (q) : 

q 



g2 



i^2{q,Po) +T^i{q.Po) 



(49) 

(50) 
(51) 
(52) 

q- q^ 

The term ciq'^ comes from the curvature of the uncoupled dispersion surface for the upper 
channel. The term is the leading order effect of the coupling. It enters the WKB 

solution at order e|r/p. The terms with negative powers of q only become large for small q, 
but they cause the WKB phase to diverge in the mode conversion region: 

(53) 



~I4 

v\ 



ciq + |?7| (a2 - bi) H ^(ai - 62) + 



C2\V\^ 



ipiiq) = Ai{q) exp J po{q') + epi{q') dq 
= Ai{q) exp |— z|f/|^ ln(g) 



Ciq-" 



+ |?7| (02 - hi)q 



I ~I4 



ai 



C2\v\^ 



(54) 



q 

The WKB amplitude is computed from the derivative of the eigenvalue Aq which 
approaches Dn for |g| ^ 1, and which is the generator for ray evolution in the upper channel: 



A. 



D 



11 



D 



22 



sgn(g)^ 



22 



D 



11 



(55) 



The p derivative of this has order e terms because of the order e terms in Dn and D22, and 
also because of the order e term in p{q): 

1 sgn(g) / q + p 



Aq |p( 



1) 



v/(p + g)2 + 4|r7 



n|2 



sgn(g) 



po(g)+epi(g) 

{q + p){dpV2-dpVi) 



2\/(g + p) +4|r/ 



7^\2 



{q + pY(V2-Vi 
{q + p)^ + 4|f/P 



(56) 



Po(i?) 



As a first step in evaluating this derivative, we can find the leading order effect in e of the 
Po{q) + ^Pi{q) evaluation: 

sgn(g) / q + p 

Po{q)+tpi{q) 



A/(j9 + g)2 + 4|?7|2 
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1 
2 



sgn(g) 



q + Po + epi 



(57) 



This gives the expansion of dpXa in e. Further simplification can be achieved by considering 
the denominators: 



[(g + Po)' + 4|r^p]-/2 



I ~I2\ 2 

\v\ 



4|r/|2 



-k/2 



q' - 2q 



~I2\ 2 

P7r 



-k/2 



q + 



I ~I2\ 2' 



-1 ~k/2 



(58) 
(59) 
(60) 



I ~I4 

\ri\ 



2 



K 



\m 

q6 



+ ... 



(61) 



Since we are working with the WKB approximation, q ^ and we can truncate this 
expansion. The order at which it can be truncated will be determined by the highest power of 
q by which this term is multiplied. For example, when multiplied by piq^, which is 0{q^), we 
need to keep terms through 0{q^^). 

Inserting all of this, including the expressions for pi, Vi, and V2, back into the derivative, 
and gathering terms of like powers in q gives 



dpXa 



q2 



+ ebiq + e|f/p (62 - h - la-i - 2ai) q^^ + C(e^ q 



(62) 



The first three terms are the most significant in the limit we are considering. So, the WKB 
amplitude can be written 



M{q) 



dK{q,p) 


-1/2 


dp 


p=p(q) 



-1 



|?7p 

g2 



+ e&ig 



'1/2 



(63) 



The calculation for the coupled WKB mode in the lower channel proceeds along similar 
lines, and uses the p representation of the equations. The result is 



V'2(p) = M{p) exp \ i\r]\^ ln(p) + ie 
with 



~|4 



ri\ (ai - h2)p (02 - hi 



p 



Ci\fj\^ 

3p3 



,(64) 



A2{P) 



?7p 



eh2P 5- + . . . 



-1/2 



(65) 



We will use these expressions for the coupled WKB solutions when matching the local 
solution to the propagating modes. 
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4.2. Local Coupled Solutions 

Now that we know the form of the incoming and outgoing WKB waves, we need to solve 
the system of equations locally in order to find the scattering coefficients which connect the 
incoming to outgoing waves. We will find the higher order corrections by expanding the local 
fields in e, and then using the 0{e) equations to get the local solution. Write the fields as 
asymptotic expansions in e: 

^i(g) = g-l'^P(l + eei(g) + 0(e^)), (66) 
Mq) = + e02(g) + 0{e^)). (67) 

Because of the form of the equations, it is convenient to expand 6i and 62 as power series in 
q: 

00 00 

n=— 00 n=— 00 

The coefficients s„ and s„ are constants which depend on e. 

In order to compact the notation, and simplify the calculations, define as the 
coefficient obtained when acting with p on one of the terms in the series expansion above: 

= i(2|f/|2-n)g-'l^l'+"-^ (70) 
= (71) 

Our equations involve the operators p, p^,q^, and Using the definition of above, 
we can evaluate the action of these operators for an arbitrary term in our series: 

Pq-Aft?+^ = a^q-Aft?+'--\ (72) 
Applying this formula twice gives us 

= (73) 

Similarly we can calculate the action of the symmetrized operator pg: 

fqq-M^+r^ = ]^{pq + qp)q-'\^\^+'' (74) 

= ^(2gj5-i)g-'l^l'+" (75) 

- (76) 
We will also use 

X By we mean the operator whose Weyl symbol is pq. Since p and q commute as variables in phase space, 
we could also have written this operator as qp 
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We can now write out our system of equations in (|39l ) using the coefficients a„ and our 
series expansions for the fields. Keeping only terms of order e, we get 

^iiV'i(g) + #2(g) = + - Ir/I' = (78) 

n n 

and 

r/^V'ilg) + ^22^2(g) = ^n^""'''^" - ^^2(g,p)g-'i'^i'-' - Sng-*i'^i'-'+" = o. (79) 

n n 

Now evaluate the operators in these equations using the expressions derived above: 

n 

+ai (ao«_i) q-'\^\'-' + h («o - ^ q"'^^^" + Cig-^l^l'+^ = 0, (80) 



-i\fi\^+n 



n 



-C2 (a-ia_2) g-'I'^l'-^ - 62 - ^ j g"*''^!'-' - a2g-*l'^l'+^ = 0. (81) 

These equations can now be solved for the coefficients s„ and s„, by matching like powers 
of q. For n ^ {—3, —1, 1, 3} these give s„ = s„ = 0. Otherwise, we have a set of linear 
equations which can be solved (using the definition of a„ to simplify the coefficients) to give 

(82) 
(83) 
(84) 
(85) 

&2 

(86) 
(87) 
(88) 
(89) 

These coefficients give the corrected local solutions in the q representation, which can be 
matched to the WKB solutions for the upper channel. For the lower channel, we need the local 
solution in the p representation. The calculation for the transformation to the p representation 
involves integrating functions which have singularities at the origin, just as we had when 
calculating the conversion coefficient in Equation (l32l) : 

ij2{p) = I dqe^'P'^Mq) (90) 



^27rz 



S3 


ICi 

~ T 




S3 


ici 




Si 


= i{bi - 02)^0 + 


bi 
2 


Si 


= i{bi - 02)0:1 - 


bi 
2 


S-1 


= i{b2 — ai)aoa^ 


1 + 


5_i 


= i(b2 — ai)aoO_ 


1 — 


S-3 


iC2 

= — aoa-i"-2 




5-3 


iC2 

= — a_ia_20_'? 
3 
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= ^L= I dqe-'''\-frq-'\'<^"~^){l + eQ^i^q) + 0{e^)) (91) 
V27r2 J 

= ^L= ! dqe-'P%-fi*q-'\^\^-^){l^ey"~Snq'' + 0{e^)). (92) 

This integral can be evaluated (by putting it into the form of Hankel's contour integral for the 
complex Gamma function), and the answer written as the power series 

Mv) = -P^'^'' f 1 + e E ^"P" + ^(^') ) ' (93) 



where terms in the series can be shown to be 

ao«i ■ ■ ■ a-n-is-n, n <0 

0, n = (94) 
'-^ , n>0. 

Inserting the values of s„ from above gives 

o-_3 = — a2«iao (95) 

5-_i = i{bi - a2)aiao - (96) 

h 

ai = i{h2 - ai)ao - — (97) 

These coefficients, together with the coefficients s„ above, give the new local solutions, which 
include the effects of the quadratic terms. 



4.3. Matching the Local and Far-field Solutions 

4.3.1. Matching the Upper Channel In order to find the connection between the incoming 
and outgoing wave fields, we need to match the WKB solutions in Section |4~n with the local 
solutions in Section 14.21 This matching will allow us to find the scattering matrix for this 
mode conversion, which determines the amplitudes and phases of the outgoing waves given 
the amplitudes and phases of the incoming waves. The various WKB waves can be identified 
with the dispersion surfaces from which their phases are calculated. This lets us label the 
branches of the dispersion surface by the type of WKB mode which it generates, as in Figure 

m 

Matching a WKB wave incoming in the upper channel requires comparison of the WKB 
and local solutions for the upper channel. The WKB solution is given in Equation (|54|) and 
the local solution is written as the series in Equation (l66l) : 

and 
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Here, the terms with negative powers of q have been dropped, since the effect of these terms 
is localized to the mode conversion region. 

The incoming matching point qMi < is now chosen where both of these expressions 
are valid, and where the neglected inverse powers of q are actually negligible. The amplitude 
and phase of the incoming WKB mode at this point are used to set the incoming amplitude 
and phase of the local solution; 

, (WKB) / X 

v-;(.)-^ ^o...,,'" ^'i"""(^) (101) 

i^i mil) 

= (102) 

Here A is the complex amplitude of the incoming wave. As in the case of the linear order 
problem in Section[3l we must choose a branch cut for the term e^*''^'^ This choice will let 
us write the matched local solution in terms of the magnitude of q. 



( 1 + ie\fi\\a2 - h)q + ^ + + ^(e^)^ for g > 



2 3 

Now choose a matching point q^o > where this local solution can be matched to the 
outgoing WKB wave. We then can write the outgoing field as 

[qMo) 

= A\qr^^^' (l + ^ema, - h)q + + + oie\ q^')^ (105) 

This expression shows that the transmission coefficient r remains the same as it was in the 
case of the linearized problem described in Section [3l The quadratic order terms introduced 
corrections to the local solution which makes the local solution better match the WKB 
solutions. 



4.3.2. Matching the Lower Channel Use the p representation to match the lower channel. 
Match the WKB solution from Equation (|64l ) , 

^^^''\P) = P'^'^' (l + ^ema, - b,)p _ f^! + !!|^ + (106) 
to the local solution from the series given in Equation (|93] ). 

.(Local)/ X P* i\fi\2 f . I ^,2/ , N ^hP , ieC2P^ . 2 _iA /\c\n\ 

'{p) = —p™ \l+te\r]\{ai-b2)p — + ^— + 0{e,p )j, (107) 

where, as for the upper channel, we have dropped the terms which involve inverse powers of 
p, since they are negligible outside of the mode conversion region. The matching proceeds 
as in the first channel, except that the initial amplitude and phase are given by the incoming 
data in the upper channel. Since the local solutions given in Equations (l66l) and (l67l) have the 
correct relative amplitude and phase, we can use the overall initial amplitude and phase from 
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( 11031 ) as the initial amplitude and phase for the lower channel. This gives us the local field for 
the lower channel; 

ij'M = ATijf°''''\p). (108) 

We now pick a point Pa/ > where we match the local solution onto the outgoing WKB 
wave. The outgoing wave is then 

, (Local) / \ 

= Ar!t2^f^i,r^\p) (109) 

W2 VPm) 



A(3*p 



(i + u\vna, - _ + !!|p! . (1,0) 

This shows that the conversion coefficient is unchanged, and that the local solution now 
matches the WKB solution much better. 



5. Comparison with Numerical Solution 

In order to demonstrate that the corrected solutions derived above do actually match better, 
we compared both the corrected local fields and the matched WKB waves to the results of 
numerical simulations. In order to avoid the singularities in the solutions, the numerical 
calculation was carried out in the Q representation, where the phase space coordinates (Q, P) 
are those defined by the linear canonical transformation in Equation (Hj). The associated 
metaplectic transformations of the coupled WKB waves from Equations (|54|) and (|64|) are 
[without the divergent inverse powers in the phase] 

4™)(g)=/ e^^^^^'^^^-=L====^dg, |g|>l (111) 

and 

^(™)(g)=/ e^^^(«'-)^^^_==^dp, H>1, (112) 

where Fi{Q,q) and F2{Q,p) are the generating functions for the linear canonical 
transformation: 

FiiQ, q) = \{Q^- 2V2Qg + g'), (113) 
^2(g,p) = -^(Q'-2v^Qp + p2). (114) 
The conjugate variables are given by derivatives of the generating functions: 

In order to evaluate the metaplectic integrals in equations (II 111) and (II 121) . we write q in 
terms of Q and possibly also derivatives with respect to Q, i.e. we find the Q representation 
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of the operator q. The properties of the generating functions allow us to do this. First, use 
Equations ^ and (II 151) to write 

g = -^(g - P) = -^(g - idQF,{Q, q)). (117) 

Now, notice that we can combine the derivative dqFi^Q, q) with the phase in the integral to 
write 

5Qi^i(g,g)e*^'^^''^ = -i9Qe*^^('3,<?)_ (118) 

We can now consider the higher order terms in the integral as a psudeodifferential operator 
acting on the function by making the substitution 

q-^l^{Q + idQ), (119) 

wherever q appears in the higher order terms. In order to do this, write the higher order terms 
as a series: 

^^^^^ ^ /^ h 1-12/2 = (120) 

oo 

n=— oo 
oo 

^ E ^(Q + '^qT (122) 

n=— oo * 

= 5i fi=(g + z9Q)^ (123) 



.V2 

We can now use this to find the Q representation of the field: 



/oo 
e^F^(Q^'^)S,{q)q-^\^\'dq (124) 
-oo 

= ^1 (^-^(g + iOq)^ P e'^^(Q,i)q-i\v? dq. (125) 



The integral here is the same as the integral in Equation (iTSl) . and evaluation gives a parabolic 
cylinder function. However, we are interested in the WKB limit, where Igl 1 and \Q\ ^ 1. 
So, we can evaluate the integral using the stationary phase approximation, which gives: 

e^F,iQ,g)^-i\n\^ ^ yi|g^V4e-.QV2-.|rypln{v^Q) ^ 0(g-2) ^ ^^^iQy2^i\rj\^ H^%^ 

Since the WKB approximation only valid far from the mode conversion, we will drop 
all terms in Si{q) with negative powers of q. Also, since the quadratic term in the phase 
varies more rapidly than the logarithmic term (for the values of Q we are considering), we 
will consider the action of the derivatives only on the quadratic term. This gives us 

V;i™)(g) ^ATs/ -^(Q + idQ)] e-QV2-.|.;Pin(V2Q) (127) 



_^g-^|riPln{v^Q) |^_l^(g ^ ^^^)^ ^-iQV2 (j^g) 



Quadratic corrections to the metaplectic formulation of resonant mode conversion 24 



gjeci(V2Q)3/3+je|^|2(a2-6i)v^(3 



^l + e(-6i + 2ci)y2g 



The appearance of the Ci term in the amplitude may be unexpected. However, under the 
Fourier transform, the pure phase exp(zcig'^/3) turns into an Airy function, which has 
amplitude variations. Converting from the q representation to the Q representation is done 
with a metaplectic transformation, which is a sort of "partial" Fourier transform. Therefore, it 
should not be too surprising that the cubic phase in q would give rise to an amplitude variation 
in g. 

A similar analysis can be computed for the lower channel. Since our corrections are 
written in the p representation, we need to make the substitution 

p=±{Q + P)^±{Q-^^Q). (132) 
Therefore, the WKB mode in the lower channel is 

= AA^e^l'^l^'-^^^Q) (^1 + e (-| + m^a^ - h)^ {V2Q) + ^ ((v^g)^ - 3zv^g)^ eff^^3 



gieC2(V2Q)V3+je|^|2(ai-fe2)V2Q 



,J 1 + e{b2 - 2c2)V2Q 

We can now compare our analytical expressions (Equations ( 11311 ) and (11361 )) and 
numerical simulations. As seen in Figure (HJ, these corrected solutions correspond closely to 
the numerical simulations. The amplitude of the analytical solutions now contain the square 
root variation which is due to action conservation, so they match the WKB solutions over a 
much wider range than before, cf. Figure ([3]). The phases of the solutions also show good 
agreement with the numerical simulations. Notice that the rapid quadratic variation in the 
phases as been subtracted from these plots, in order to better show the effects of the coupling 
and the higher order terms. 



6. Conclusions 



In this paper we have shown how to extend the metaplectic formulation of resonant mode 
conversion to include the effects of quadratic order variations in the dispersion matrix. A 
corrected local solution was derived, and matched onto far- field WKB solutions, showing that 
the transmission and conversion coefficients are unchanged by the new quadratic order terms. 
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The corrected solution also matches the far-field solutions over a much larger region, as is 
illustrated by numerical simulations. 

Appendix A. The Normal Form in One Spatial Dimension 

These three appendices will demonstrate that, for one-dimensional problems, the dispersion 
matrix can be put into normal form order by order in powers of the phase space variables. In 
the present one-dimensional setting, by 'normal form' we mean that the off-diagonal terms 
are constant and that at linear order the new diagonal terms still form a conjugate pair. In 
two, or more, spatial dimensions the normal form is characterized by the requirement that the 
diagonals Poisson-commute with the off-diagonals, which ensures that the off-diagonals are 
constant following rays generated by the diagonals. This will be discussed elsewhere. 

In order to put the matrix into normal form, polarization vectors must be chosen so 
that the off-diagonal elements at higher order vanish. The strategy for this calculation is 
broken into three parts for clarity. The first step is to explicitly calculate the transformation 



which puts the matrix into normal form through second order (Appendix A I as a way to 
introduce key ideas in the simplest setting. The second step is to outline the calculation which 
would be needed to put the matrix in normal form to arbitrary order, and show that there are 



sufficient free parameters in the transformation to achieve normal form (Appendix B I. In both 



Appendix A and Appendix B we ignore Moyal corrections, which significantly complicate 



the picture. The final step is to examine the effect of the Moyal corrections (needed when 



multiplying symbols of operators) which were neglected in the first two steps (Appendix C). 

As described in Section[3l and shown in lfT8l [6l. the 2x2 symbol of the dispersion matrix 
can be put into the following "normal form" at linear order: 

DNF(g,p)= ( Tf ). (A.l) 

y T] q J 

Here, is a constant since we are working in a two-dimensional phase space. The higher 
order corrections to this matrix appear at quadratic order in the phase space variables: 

D(g,p) = T>Mq,p) + e'D2(g,p) + 0{e'), (A.2) 

where e is a formal parameter introduced to keep track of the ordering. Each element of D2 
can contain terms which are quadratic in the phase space variables z = {q,p). When needed, 
we will also write the first order and second order terms by displaying the monomials: 

Di (g, p) = gDg + pBp, D2 (g, p) = g^Dg^ + pgBpg + p^Dpp, (A.3) 

where each of these constant 2x2 matrices of coefficients is hermitian and, therefore, 
generically have four real parameters. The first order matrices are particularly simple, but 
the general second order terms have twelve real parameters. 

Because D is in normal form to linear order already, we will use a near-identity change 
of polarization basis to bring it into normal form at second order. We write Q as 

Q(g,p) = 1 + eQi = 1 + epQp + egQg, (A.4) 
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where Qj, and are constant 2x2 complex matrices. Hence, we have eight complex 
parameters to work with, or, equivalently, sixteen real ones. Carrying out the congruence 
Q^^DQ, collecting orders in e, we find: 

D' = Do + e [Di + QIDo + DoQi] + + QIDi + DiQi + QIDoQi] + 0{€^){A.5) 

The sixteen real parameters in Q must be chosen so that the off-diagonals of the bracketed 
terms are zero. The term Di is already diagonal, so at (9(e) we have the two conditions: 

[QtDo + DoQ,] ^2 = 0, [QjDo + DoQj = 0. (A.6) 

(The other off-diagonal terms are the complex conjugates of these expressions because 
the congruence preserves the hermiticity.) Each of these expressions consists of two real 
conditions on the sixteen parameters in Qi, hence four of our degrees of freedom are used up. 
These conditions are particularly simple because Dq is so simple: [Q^]^]^ + [Qq]22 ~ ^ ^^'^ 
[Qplii + [Qp] 22 ^ ^- ^(^^) have the three conditions: 

[D,, + QjD, + D,Q, + QjDoQ,] = 0, (A.7) 

[D,p + QjDp + DpQ, + QtD, + D,Qp + QjDoQp + Qj,DoQ,] = 0, (A.8) 

[Dpp + QjDp + DpQp + QjDoQp] = 0, (A.9) 

which use up another six parameters. Thus, to ensure that that the off-diagonal terms are 
constant to order e'^ we use ten of the sixteen parameters in Qi. Notice the pattern: at order 
e"* there are m + 1 monomials (g*", q"^~^p, . . . p^). Hence, at each order we require m + 1 
complex (2m + 2 real) conditions to be satisfied. Normalizing the Poisson bracket of the 
diagonals to linear order so they form a conjugate pair takes one more parameter. 

In summary, by counting parameters, it would appear we can generically put an arbitrary 
2x2 dispersion matrix that is quadratic in the phase space variables into normal form. We 
now perform an explicit calculation choosing a particular parameterization to demonstrate this 
concretely. To simplify notation, we write the off diagonal elements of the second order terms 
as [Dpp]i2 = dpp, [Dqp]i2 = dqp, and [Dqg]i2 = dqq. Some algebra shows that the matrices Qg 
and Qp can be parameterized in the following manner: 



Q,= l "V-^-^" + ^-) ^- + "* |. (A.10) 

-dlq -e-'^^-dqq{a* + dqp) 

Qp= I 1, (A.ii) 



a 



-e-''^^^^p 



Here, is the phase of the coupling, f] = \ fj\e'"^, and a is given by the solutions to the equation 

a*^ + dgpa* + dppdqq = 0. (A. 12) 

This transformation will make the off-diagonals in D' constants plus terms starting at 0{e^). 

In order to obtain the normal form of D, we need to normalize the first order terms in the 
diagonals to have unit Poisson bracket. We can interpret the new expressions at linear order 
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as new canonical variables {q',p') if they are related to the old {q,p) by a linear canonical 
transformation: 





(A.13) 



~ |i3|V2 23^(^3^/*) 1 + 2^{afi* 

Here the normalization is introduced to make {q',p') a canonical pair. It is given by the 
Poisson bracket of the transformed diagonal elements, 

B = {D[„ D',,} = (1 - 4(3ft(ar/*))2 - A^idiV^Mdsfj*)) {q,p} (A.14) 
= {q,p}{l + 0{e')). (A.15) 

Appendix B. Extension to Arbitrary Order 



The calculation of Appendix A can be extended to put the dispersion matrix into normal 
form order by order in powers of the phase space variables. In this appendix, we outline 
the transformation, and show that there are enough free parameters in the near-identity 
transformations at each order to eliminate all non-constant terms in the off-diagonals elements 
of the dispersion matrix. We still ignore Moyal corrections in this calculation. Those will be 
discussed in the next Appendix. 

Start with the a dispersion matrix that is in normal form through order A^. This means 
that we can write the matrix as 

D(g,p) = e°Do(g,p) + e^Di(g,p) + e^B2{q,p) + e'D3(g,p) + . . . 

+ e^Djv(g,p) + e^+iD^+i(g,p) + . . . , (B.l) 
where Dq is the constant coupling 

Do(g,p) = ° ly (B.2) 

and each of the matrices T)j{q,p), for j < N, are diagonal matrices with entries which are 
homogeneous polynomials of order j. At 0{e'^) there are m + 1 such matrices, one for each 
monomial for m' = 0,1, . . .m. 

In particular, we can assume that the first order term has been transformed into the normal 

form 



Di(g,p)= / . (B.3) 




We now want to apply a near-identity change of polarization basis which puts this matrix 
into normal form through order + 1. Write the transformation matrix as 

Q(g,p) = l + e^Q(g,p), (B.4) 

where Q(g,p) is a matrix of homogeneous polynomials of order A^. There are A^ + 1 
monomials of this order, hence there are 4(A^ + 1) complex parameters (8(A^ + 1) real 
parameters) to work with. 
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The transformed dispersion matrix can then be written as 

D' = D + e^q} ■ D + e^D ■ Q + e^^Q^ ■ D ■ Q. (B.5) 

The last term in this expression starts at order 2N . Since we only need to consider terms of 
order and + 1, we can drop the last term except in the case = 1. However, this is the 
case considered in Appendix A, and so we do not need to consider it here. We can now use 
Equation (|B.1I) to group the remaining terms in D' by their order; 

D' = e°Do + e^Di + . . . + e^"^D^„i C(< A^) (B.6) 

+ e^(D;v + ■ Do + Do ■ Q) 0{N) (B.7) 

+ e^+^(D^+i + Qt ■ Di + Di ■ Q)0(Ar + 1) (B.8) 

+ ... C(>A^ + 1). (B.9) 

The order A^ and A^ + 1 terms are the ones which now need to put into normal form. For this 
1 -dimensional problem, this means that the matrices 

Q^Dq + DoQ (B.IO) 

and 

D^+i + Qt-Di + DiQ (B.ll) 

must be put into diagonal form. Using Equations (IB. 21) and (IB. 31) . we can simplify the 
equations for the elements of Q. The constraint on the order A^ matrix gives 

[Q]ii + [QI22 = 0. (B.12) 

These are A^ + 1 complex (2A^ + 2 real) conditions, one for each monomial of order A^. This 
leaves us with (8A^ + 8) - (2A^ + 2) = 6A^ + 6 free parameters. 
The C(e^+^) normal form condition gives 

p[Q]i2 - g[Ql2i = [D^+i]i2, (B.13) 

assuming that D^r+i is a hermitian matrix. There are A^ + 2 monomials of order A^ + 1. 
Therefore, this set of constraints uses up only another 2A^ + 2 real parameters out of our 
remaining allotment of 6A^ + 6. Therefore, generically, there is more than enough freedom to 
carry out the normal form transformation, order by order. We now consider how this picture 
changes when we include Moyal corrections. 



Appendix C. Moyal Corrections for Phase Space Dependent Changes of Polarization 

In the previous section, a phase space dependent change of polarization is used to put the 
dispersion matrix into the form where its off diagonal elements are constants with a small 
perturbation that starts at order e"^. This transformation is achieved through the conjugation 

D'(z) = Qt(z) ■D(z) ■ Q(z). (C.l) 

However, since the matrices D and Q are actually matrix valued symbols of operators, we 
really need to use the Moyal star product to multiply the elements of the matrix. The 



Quadratic corrections to the metaplectic formulation of resonant mode conversion 



29 




noncommutative star product is used in the symbol calculus so that the symbols (functions 
on phase space) maintain the commutation relations of the original operators. So, we should 
actually use the expression 

D'(^) = Qt(^)*D(2)*Q(^), (C.2) 

where the matrices are multiplied in the usual way, but the elements of the matrices are 
multiplied using the star product. In this section we will first consider how the Moyal 
corrections mix orders in e. We will then use these results to sketch an argument that shows 
we can carry out the normal form transformation order by order, including Moyal terms. All 
infinite series expressions should be interpreted formally. We do not consider convergence, 
nor do we worry about the asymptotic character of these expressions. 

The Moyal star product of two symbols A{z) and B{z) is often written as a formal power 
series (see |[T9ll with h set to 1 since we are studying classical fields). In general, if the symbols 
are transcendental functions, this series will contain infinitely many terms: 

A{z)*B{z) = A{z)exp | - — Ja/s— ] B{z) (C.3) 

B{z) (C.4) 

= A{z)B{z) + Ua,B]-\ {dlAdlB -2dgdpAd,dpB + dlAd^B) + .(C.5) 

Here, Jap is the symplectic matrix. In the last line, the conjugate phase space coordinates 
z = (q, p) are written explicitly to illustrate the nature of the terms in the previous sum. 

We are assuming that the symbol of the dispersion matrix has a well defined Taylor's 
series in the mode conversion region, and that the lowest order terms in the series dominate. 
This can be expressed mathematically by introducing the small parameter e. Using the multi- 
index notation, we can expand the symbols in terms of all possible mononials: 

oo oo 

A{z) = Y,e'' Yl ^rnz""^ ^(^) = E E ^^-^^ 

M=0 \m\=M M=0 \m\=M 

The multi-index m is a pair of integers (mi, 1712) with \m\ = mi + m2, and z"' defined as 

The terms amz"^ and bnz"^ are homogenous polynomials of order m in z. We will consider 

their dependence on e momentarily. The star product of the series A and B is: 

00 

A{z)*B{z)= E 5^«-^n2™*^". (C.8) 

N,M=0 \m\=M\n\=N 

Now consider the star product of two generic monomials, z"* * z". Because the star product 
can be written as a series in powers of derivatives, the star product 2;™ * 2" contains powers of 
z ranging from \l\max = m + n to \l\min = \{m + n) — 2min{m, n)\. All of the coefficients 
of this polynomial can be calculated from equation (|C.3I) : 



+ Y ci{m,n)zK (C.9) 
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We will not need their explicit form, simply the fact that they are well defined functions of 
(m, n) and /. It is also important to note that, because each term in the Moyal series (IC.3I) 
involves derivatives acting both to the left and right, the series (IC.8I) descends in steps of order 
2 in 2 (e.g. |m|, |m| — 2, |m| — 4, . . .). This result is significant for us because it means that the 
star product only introduces monomials of degree less than the degree of the ordinary product. 
We can use this result to write 

oo L— 1 

A{z) * B{z) = A{z)B{z) + E E E (C.IO) 

N,M=0 \m\=M \n\=N \l\=0 

oo L— 1 

= A{z)B{z) + Y,^'^J2^'i'' (C.ll) 

L=0 \l\=0 

(oo L-1 \ 

E E (^^n • (C.12) 
L=0 |«|=0 / 

This means that the original functional form we assumed, which ties the powers in e directly to 
the order must be modified. We replace it by the assumption that the coefficient of e™ has 
^-dependent coefficients that include no terms of higher order than \m\, but which can include 
all lower order terms in z. It is clear that the Moyal corrections significantly complicate the 
algebra we have to deal with. In spite of this, it is still possible to show that we can arrive at 
well-posed formal iterations schemes for constructing these series. 

As an example, consider the following problem: suppose we are given the symbol 
B{z) = 1 + Em=i where the bm (z) are independent of e and include terms in z 

of maximal order |m|. We now ask if we can find a symbol A(z) of similar form such that 
A* B = 1. Write the formal series A{z) as a Moyal product of terms: 

A{z) = . . . (1 + e'^asiz)) * (l + e^a2{z)) * (1 + ea,{z)) . (C.13) 
Acting with (1 + eai(2;)) we have: 

1 = (1 + eai{z)) * (1 + €bi{z) + e%{z) + e%{z) . . .) . (C.14) 
We collect terms in e as usual, but now we must respect the Moyal product ordering: 

1 = 1 + e [ai{z) + bi{z)] + ^z) + ai{z) * bi{z)] + .... (C.15) 

This fixes ai(z) = —bi{z) and determines all of the higher order terms, which now become 

b'^{z) = bn{z) — bi{z) * 6„,„i(z). At next order, acting with (1 + e'^a2{z)), we have: 

1 = (1 + e^a2{z)) * (1 + eX{z) + ...)• (C.16) 

This fixes a2{z) = —62 (-2) and modifies all of the higher order even powers in e. Clearly this 
procedure can be carried out to arbitrary order. Using the insights gained from this simple 
example, we now return to the problem of casting D into normal form. 

Once again, we assume that 0(2) is in normal form up to 0{e^). To be precise: we 
assume that 

00 

D(^) = E E ^"'^n^{z), (C.17) 
M=0 \m\=M 
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where the terms up to are in normal form. The 2x2 hermitian matrices Y)ra{z) are 
now allowed to contain arbitrary terms in z up to at most We use a congruence 

transformation Q of the same form as (IB.4I) . but we must now allow Q to be a general series 
in z of terms up to order |m| = A^. The calculation proceeds much as before, but now the 
Moyal product is used instead of the ordinary product when two matrix entries are multiplied. 
The terms up to order — 1 are unaffected. At 0{e^) we have 

D'^(^) = D^(2) + Q^(^)Do + DoQ;v(^). (C.18) 

Because Do has no 2-dependence, this condition is unchanged from the case where Moyal 
corrections are ignored: 

[Qiv]llW + [Q;r]22(^)=0. (C.19) 

At 0{e^^^) we find now two conditions: 

[Q^]21 * g - P * [QAr]l2 = -[D,v+l]l2, (C.20) 

and 

q * [Q^]i2 - [Q^]2i * P = -[D,v+i]i2, (C.21) 

where we have used [DAr+i]i2 = [D7v+i]2i by hermiticity. In addition, we require that the 
new diagonals at linear order in q and p are a conjugate pair, which is one further condition. 
At a fixed order N in z, there are + 1 monomials. Hence, because these expressions include 
terms from z^ up to order A^, there are 

A/=0 

conditions that must be satisfied in each of these three equations (IC.19I) . (IC.20I) and (IC.21I) . 
This makes for 3(A^^ + 3A^ + 3) real conditions in total at order A^. But, there are 
4(A^2+3A^+3) parameters to work with (the four entries of Qn{z) each have (A^^+3A^+3) /2 
monomials with complex coefficients). Therefore, formally, the normal form conditions can 
be satisfied order by order in e and z through a proper choice of a z-dependent choice of 
polarization basis, including Moyal corrections. 
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